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AEROELASTIC  RECIPROCITY:  AN  INDICIAL 
RESPONSE  FORMULATION 


Introduction 

The  present  report  describes  a  continuation  of  research  initiated  and  reported  under  grant  AFOSR-89-0502. 
The  focus  of  the  research  has  been  an  experimental  and  analytical  investigation  of  airfoil  indicial  response 
aerodynamics  as  related  to  the  convolution  integral  formulation  for  unsteady  airfoil  loading. 

The  convolution  integral  formulation  for  the  lift  on  a  flat  plate  airfoil  in  arbitrary  motion  has  been  derived 
in  Reference[l].  The  formulation  requires  an  indicial  lift  function  which  by  definition  is  the  transient  lift 
response  to  a  step  change  in  angle  of  attack.  The  indicial  lift  response  for  a  flat  plate  airfoil  experiencing  a 
step  change  in  angle  of  attack  due  to  plunge  was  first  derived  by  Wagner.  In  Reference  [1],  Wagner's  function 
is  derived  under  the  assumption  of  small  perturbations  using  a  Fourier  integral  of  Theodorsen's  function  for 
harmonic  motions.  At  the  onset  of  the  step  Wagner's  function  jumps  instantaneously  to  one-half  the  steady 
state  lift.  Also  derived  in  Referencefl]  is  Kussner’s  function  for  the  lift  on  a  flat  plate  airfoil  entering  a  sharp- 
edged  gust.  The  initial  lift  in  this  case  is  zero  and  approaches  steady  state  as  the  gust  overtakes  the  airfoil.  In 
the  case  of  a  2D  airfoil  of  finite  thickness.  Reference  [2]  shows  that  for  an  airfoil  having  a  trailing  edge  with  a 
finite  angle  of  convergence  the  indicial  lift  at  the  step  onset  is  zero.  This  suggests  that  the  subsequent  transient 
lift  for  finite  airfoils  at  small  incidence  may  be  similar  in  form  to  Kussner’s  function.  Reference  [2]  cites  other 
numerical  studies  for  finite  airfoils  which  appear  to  agree  with  the  zero  initial  lift  result  Reference  [3]  uses  the 
Laplace  transform  method  to  derive  the  lift  transfer  function  for  a  number  of  airfoil  motions.  The  Laplace 
domain  formulation  includes  the  indicial  lift  function  and  results  are  presented  using  Kussner's  function  to 
describe  the  indicial  lift. 

Reference  [4]  describes  a  tow  tank  study  designed  to  measure  the  normal  force  response  of  a  NACA  0015 
airfoil  experiencing  sudden  step-like  changes  in  angle  of  attack  by  rotation  about  the  quarter.  The  motivation 
for  these  experiments  was  to  study  airfoil  indicial  response  aerodynamics  as  defined  in  the  theory  of  nonlinear 
mathematical  modeling  for  aerodynamic  systems  [5].  These  experiments  involved  strain  gauge  load  cell 
measurements  of  the  transient  normal  force  loading  on  an  airfoil  undergoing  a  sudden  change  in  angle  of  attack 
of  approximately  Aa  =  ±1°.  The  angle  of  attack  prior  to  the  step  onset  (a^)  was  steady  and  was  varied  from 
run-to-run  over  the  range  2°  <  aQ  <  60°  . 

In  the  experiments  of  Reference  [4]  the  test  rig  may  experience  large  inertial  and  aerodynamic  loading  due  to 
the  rapid  starting  and  stopping  required  to  impart  the  step.  Therefore,  an  important  issue  is  the  degree  to  which 
aeroelastic  reactions  deform  the  structure  thereby  influencing  the  output  of  the  strain  gauge  bridge.  Knowledge 
of  these  reactions  is  useful  in  comparing  these  strain  gauge  data  with  the  theoretical  indicial  responses  of 
Wagner  and  Kussner.  Part  I  of  this  report  describes  an  aeroelastic  analysis  of  the  Ohio  University  tow  tank 
test  rig.  The  model  is  based  on  the  mode  superposition  method  for  structural  systems  and  classical  linear 
airfoil  theory.  The  Laplace  transform  method  is  used  to  solve  the  equations  of  aeroelasticity  in  closed  form. 

Part  II  of  this  report  describes  an  inverse  aeroelastic  analysis  designed  to  determine  the  aerodynamic  indicial 
response  from  the  aeroelastic  response.  The  method  provides  a  theoretically  sound  procedure  for  determining 
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the  indicial  response  experimentally,  even  though  the  indicial  response  is  defined  mathematically  as  the  response 
to  a  step  motion  of  diminishing  amplitude.  Part  III  describes  indicial  responses  computed  with  the  method  of 
Part  II  using  tow  tank  strain  gauge  data. 

Ohio  University  Tow  Tank 

A  schematic  of  the  O.U.  tow  tank  is  shown  in  Figure  la.  The  facility  consists  of  a  large  tank  with  a  six 
inch  chord  NACA  0015  airfoil  suspended  vertically  in  the  water  with  a  submerged  length  of  42.0  inches.  A 
carriage  moves  in  translation  at  2  ft/s  (Re  =  9.5E4)  along  roller  bearings  fixed  to  I-beams  which  span  the  tank. 
The  airfoil  is  driven  in  rotation  by  a  drive  shaft  fixed  to  the  airfoil  quarter  chord  at  one  end,  and  coupled  to  a  3.5 
hp  stepper  motor/gear  box  apparatus  at  the  other  end.  Figure  lb  shows  details  of  the  drive  shaft  and  airfoil. 
Shown  here  are  dimensions  in  inches,  and  a  numbering  scheme  (1  through  23)  defined  for  the  purpose  of 
discretizing  the  mass  of  the  structure  as  will  be  discussed  in  detail.  The  drive  shaft  has  a  number  of  variations 
in  cross  section  over  the  length  of  the  shaft  which  must  be  considered  in  modeling  the  aeroelastic  response  of 
the  structure.  Near  the  middle  of  the  drive  shaft  is  a  machined  rectangular  section  which  has  a  strain  gauge 
load  cell  adhered  to  the  shaft.  The  load  cell  section  is  discretized  into  mass  elements  5  through  9  with  the  strain 
gauges  located  at  the  centroid  of  element  7.  The  strain  gauge  circuit  is  electrically  compensated  to  be  sensitive 
to  chord  normal  forces  only.  The  upper  most  mass  element  1  is  made  of  steel  while  all  other  parts  are 
aluminum.  The  modulus  of  elasticity  (E)  for  steel  was  taken  to  be  30x10^  psi,  and  for  the  aluminum  was 
measured  to  be  7.3x10^  psi.  The  modulus  of  rigidity  (G)  was  taken  to  be  10x10^  psi  for  the  steel  and  for 
aluminum  4x10^  psi.  The  area  moment  of  inertia  of  the  airfoil  about  the  chord  line  is  0. 12  iri^  and  the  center 
of  mass  of  the  airfoil  is  located  1.335  inches  aft  of  the  pitch  axis  at  the  quarter  chord.  The  polar  moment  of 
inertia  of  the  airfoil  about  the  pitch  axis  is  6.71  iri^.  The  densities  of  steel  and  aluminum  were  taken  to  be  0.3 
and  0. 1  lbm/in^.  The  mass  of  the  airfoil  per  unit  length  is  0. 1598  lbm/in. 


Theoretical  Linear  Normal  Force  Response  for  Small  .Angles  of  Attack 

In  an  incompressible  flow  at  small  incidence,  the  theoretical  linear  normal  force  coefficient  response  of  an 
airfoil  given  an  instantaneous  step  change  in  angle  of  attack  by  rotation  about  the  quarter  chord  is  related  to  the 
indicial  lift  function  ,  T(t),  and  is  given  by: 


C^t)  =  CNo  +  nAa  [  6(t-0)  +U,(t-0)]  +  2  n Aa  (  T(t)  -  r (t) ) 

where  C^Q  is  the  initial  normal  force,  Aa  is  the  step  amplitude,  6  is  the  time  derivative  of  the  Dirac  delta 

function  5,  andT  is  the  time  derivative  of  T.  The  first  tw7o  terms  in  brackets  in  Equation  (1)  are  generalized 
functions  [6]  which  describe  the  noncirculatory  component  of  the  loading,  while  the  last  group  of  terms  gives 
the  circulator}7  component.  For  a  flat  plate  airfoil  the  indicial  function  can  be  represented  by  a  two  pole  curve 
fit  to  Wagner's  function: 


T(t)  =  [  1  -  0.165e-*0455t  -  0.335e-3t  ] 


(lb) 


To  Pulley  on  Ceiling 

Angle  of  Attack  Nylon  Coated  Steel  Cable 
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Other  workers  have  used  Kussner's  function  for  the  lift  on  an  airfoil  penetrating  a  sharp-edged  gust,  in  which 
case: 

T(t)  =  [  1  -  0.5e-13t  -  0.5e~l  ]  (lc) 

Transient  normal  force  responses  of  a  NACA  0015  airfoil  undergoing  sudden  step-like  changes  in  angle  of 
attack  by  rotation  have  been  measured  in  the  Ohio  U.  tow  tank.  Figure  2  shows  angle  of  attack  data  for  a 
typical  run  (small  spikes  are  electrical  noise).  The  onset  angle  is  2.09  and  the  step  amplitude  is  approximately 
+ 1 .25°.  The  motion  resembles,  to  a  reasonable  approximation,  a  small  amplitude  ramp  motion.  Equation  (1) 
is  compared  below  with  the  experimental  strain  gauge  data  corresponding  to  the  motion  of  Figure  2.  To 
facilitate  the  comparison,  the  coefficient  2jc  (for  a  flat  plate  airfoil)  on  the  circulatory  part  in  Equation  (1)  has 
been  replaced  by  the  static  normal  force  curve  slope  (for  the  present  NACA  0015)  which  has  been  measured  in 
an  independent  test  [4].  This  substitution  is  necessary  for  the  response  of  Equation  (1)  to  approach  the  same 
steady  state  as  the  experimental  response. 

Part  I:  Aeroelastic  Analysis 

The  present  analysis  is  based  on  the  mode  superposition  method  for  a  forced  structural  dynamic  response 
[7],  and  linear  airfoil  theory  formulated  in  terms  of  the  convolution  integral  for  the  loading  on  an  airfoil  in 
arbitrary  motion  [1],  The  structure  to  be  modeled  has  been  shown  in  Figure  lb.  Because  the  pitch  axis 
(drive  shaft  axis)  does  not  coincide  with  the  center  of  mass  of  the  lower  part  of  the  structure  (mass  elements  13- 
23),  it  is  necessary  to  consider  the  coupling  between  the  chord  normal  and  torsional  aeroelastic  degrees  of 
freedom. 

System  Representation 

As  illustrated  in  Figure  lb,  the  structure  has  been  discretized  into  23  mass  elements  w'hich  are  considered  to 
be  concentrated  at  the  centroid  of  each  element.  Under  the  influence  of  aeroelastic  loading,  the  masses  will 
deflect  normal  to  the  airfoil  chord  referred  to  here  as  the  normal  degree  of  freedom  (NDOF),  as  w’ell  as  in  torsion 
about  the  pitch  axis  hereafter  referred  to  as  the  TDOF.  The  lowest  nine  masses  (15-23)  represent  the 
submerged  portion  of  the  airfoil,  and  masses  5  through  9  correspond  to  the  rectangular  cross  section  load  cell. 
Masses  are  concentrated  on  the  load  cell  to  obtain  good  resolution  of  the  deformation  and  strain.  In  the  TDOF, 
the  mass  elements  are  replaced  by  polar  mass  moment  of  inertia  elements.  For  the  TDOF  analysis,  the  structure 
is  discretized  in  the  same  wfay  as  the  NDOF  for  elements  1  through  12,  however,  elements  13  through  23  are 
lumped  into  a  single  inertia  element  giving  a  total  of  13  inertia  elements.  The  reduction  in  elements  in  the 
TDOF  is  based  on  the  fact  that  the  torsional  stiffness  of  both  the  mounting  block  (element  13)  and  airfoil  are 
much  larger  than  the  drive  shaft,  and  consequently  the  mounting  block  and  entire  airfoil  experience  nearly  the 
same  TDOF  deflection.  The  resulting  discretized  structure  may  be  described  mathematically  in  terms  of  a 
diagonal  mass  matrix  [M],  a  symmetric  NDOF  flexibility  matrix  [AN],  a  diagonal  polar  mass  moment  of  inertia 
matrix  [J],  and  a  symmetric  TDOF  flexibility  matrix  [AT].  The  flexibility  coefficient  AN(i,j)  is  by  definition 
the  deflection  of  mass  i  due  to  a  unit  force  applied  at  mass  j,  w^hile  AT(i,j)  gives  the  angular  rotation  of  polar 
inertia  i  due  to  a  unit  torque  applied  to  inertia  j.  The  matrix  [AN]  has  been  computed  by  assuming  the  structure 
to  deform  as  a  cantilever  beam  and  integrating  the  second  order  differential  equation  for  the  elastic  beam  curvature 
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[8].  A  condition  of  continuous  slope  is  applied  at  discontinuities  in  area  moment  of  inertia  of  the  drive  shaft. 
The  matrix  [AT]  has  been  computed  using  mechanics  for  shafts  in  torsion[9].  The  analytically  derived 
flexibility  matrices  have  been  validated  to  some  extent  by  loading  the  structure  and  measuring  the  deflection  at 
selected  points.  The  NDOF  stiffness  matrix  [KN]  and  TDOF  stiffness  matrix  [KT]  are  computed  by  inverting 
[AN]  and  [AT],  respectively,  and  has  been  done  using  an  IMSL  subroutine. 

Eigenvalue  Problem 

The  system  matrices  above  have  been  used  to  solve  the  eigenvalue  problem  for  the  natural  frequencies  and 
mode  shape  vectors  of  the  system.  The  eigenvalues  (natural  frequencies)  and  eigenvectors  (modal  vectors)  have 
been  computed  using  the  IMSL  subroutine  EVCRG.  In  the  mode  superposition  analysis,  the  two  lowest 
NDOF  modes  (NDOF1  and  NDOF2)  have  been  used,  while  for  the  TDOF  only  the  lowest  fundamental  mode 
has  been  retained.  The  corresponding  modal  vectors  are  plotted  in  Figure  3  and  given  numerically  in  Table  1. 
The  modal  vectors  are  designated  below  as  ^  f°r  NDOF1,  NDOF2,  and  TDOF, 

respectively. 

Table  1.  Modal  Vectors  for  NDOFL  NDOF2,  and  TDOF 


> 

4*ni 

4*N2 

<t>T 

1 

0.003 

-0.002 

0.0392 

2 

0.0011 

-0.0073 

0.0647 

3 

0.0025 

-0.0156 

0.0983 

4 

0.0156 

-0.0922 

0.1178 

5 

0.0229 

-0.1332 

0.2080 

6 

0.0270 

-0.1542 

0.3687 

7 

0.0318 

-0.1778 

0.5294 

8 

0.0375 

-0.2038 

0.6900 

9 

0.0439 

-0.2319 

0.8507 

10 

0.0620 

-0.3067 

0.9417 

11 

0.0928 

-0.4287 

0.9629 

12 

0.1189 

-0.5252 

0.9867 

13 

0.1416 

-0.6019 

1.0000 

14 

0.1661 

-0.6790 

1.0 

15 

0.2090 

-0.7609 

1.0 

16 

0.2803 

-0.7984 

1.0 

17 

0.3642 

-0.7455 

1.0 

18 

0.4581 

-0.6027 

1.0 

19 

0.5595 

-0.3783 

1.0 

20 

0.6661 

-0.0871 

1.0 

21 

0.7760 

0.2522 

1.0 

9 


22  0.8877  0.6202  1.0 

23  1.0000  1.0000  1.0 

*The  position  of  element  i  is  located  at  the  centroid  of  the  element  (see  Figure  lb). 

The  nondimensional  natural  frequencies  (cob/U,  b=semi  chord,  U^velocity)  for  the  NDOF  modes  were 
computed  to  be  co  =  6.374  and  0)^2  =  26.985  ■>  f°r  TDOF  to^p  =  65.782.  Notice  that  the  initial 
part  of  the  angle  of  attack  data  (the  ramp)  of  Figure  2  has  a  slope  nearly  the  same  as  a  sine  function  with  a 
nondimensional  frequency  near  50.  This  frequency  is,  in  a  sense,  the  excitation  frequency.  Because  the 
NDOF2  and  TDOF  natural  frequencies  are  near  the  excitation  frequency,  it  is  appropriate  to  consider  these  modal 
responses.  As  seen  below,  apparent  mass  and  inertia  effects  tend  to  reduce  the  aeroelastic  frequencies  even 
further. 

Structural  Dynamics  Model 

The  aeroelastic  structural  response  of  the  O.U.  test  rig  has  been  modeled  using  the  mode  superposition 
method  and  classical  linear  airfoil  theory.  For  modeling  the  rigid  body  rotation  of  the  rig,  a  rotation  degree  of 
freedom  (RDOF)  is  introduced.  The  RDOF  is  used  to  simulate  the  change  in  angle  of  attack  of  the  structure 
due  to  the  rotation  imparted  by  the  stepper  motor.  The  NDOF  and  TDOF  as  defined  above  simulate  only 
deformations  relative  to  the  instantaneous  position  of  the  top  of  the  drive  shaft  where  both  the  NDOF  and 
TDOF  deflections  are  always  zero. 

In  the  NDOF  and  TDOF,  the  deflections  of  the  structure  are  given  by  the  normalized  mode  shape  vectors 
multiplied  by  time  dependent  modal  amplitudes.  The  total  deflection  of  the  structure  is  defined  by: 

vnW  =  +  ♦N2<fc(0  ® 

vT  =  <(>Tg(t) 
a(t)  =  p(t) 

where  and  are  the  deflection  vectors  in  the  NDOF  and  TDOF,  respectively,  and  cx(t)  describes  the  RDOF 
motion.  The  quantities  ql  (t),  q2(t),  and  g(t)  are  the  modal  amplitudes  which  are  the  generalized  coordinates  of 
the  system.  The  scalar  RDOF  motion  variable  p(t)  is  the  magnitude  of  the  nominal  angle  of  attack  as  the 
structure  is  pitched  and  is  an  input  to  the  model.  This  is  the  angle  measured  by  the  rotational  potentiometer 
illustrated  in  Figures  la  and  2. 

Substituting  Equations  (2)  into  the  equations  of  motion  for  the  system  and  using  the  orthogonality  of  the 
modal  vectors  w.r.t.  the  system  mass  and  stiffness  matrices  to  uncouple  the  modal  equations  [7]  yields  the 
following  three  generalized  equations  of  motion: 

Miqi  +  Kiqj  =4>nif 
M2q2  +  K2q2  =  <t>N2^ 

Jg  +  Cyg  +  Kjg  =  <j»jT 


(3a) 

(3b) 

(3c) 
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where  the  generalized  system  properties  are  scalar  quantities  given  by: 

Ml  =  <t>N![M]<t>Nl  ,  M2  =  <|>N2lM]4,N2  ,  J  =  ‘I'tWt 

Ki  =  ♦JjIKNl+Nl  ,  K2  =  ♦n2Ekn]+N2  ,  KT  =  ^[KTI<t.T  (4) 

Cx  =  2^xJ(ox 

The  damping  term  in  the  last  of  Equations  (3)  is  included  to  account  for  a  small  amount  of  rotational  "play”  in 
the  airfoil  drive  shaft  which  tends  to  cause  oscillations  in  TDOF  to  damp  out.  In  Equation  (4),  Ly  ls 
damping  ratio  which  was  estimated  to  be  0.05  and  co  j  is  the  natural  frequency  in  the  TDOF.  The  solution  of 
Equations  (3)  for  the  generalized  coordinates  q^t),  q2(t),  and  g(t)  subject  to  a  prescribed  RDOF  input  parameter 
p(t),  is  described  below.  All  terms  in  Equations  3  have  been  nondimensionalized  based  on  freestream  velocity, 
semi  chord  length,  and  the  density  of  water. 

Aeroelastic  Normal  Forces 

The  normal  force  acting  on  the  structure  is  decomposed  into  a  rigid  body  force  vector,  FR  ,  associated  with 
the  RDOF,  an  aeroelastic  force  vector,  ,  due  to  time  dependent  aeroelastic  deflections  along  the  span  of  the 
airfoil,  and  an  inertial  force  vector,  Fj  ,  so  that  in  Equations  (3),  F  =  FR  +  FA  +  Fj  . 

The  rigid  body  force  vector  is  the  ideal  aerodynamic  loading  response  to  a  step  change  in  angle  of  attack  due 
to  rotation  about  the  quarter  chord.  These  aerodynamic  forces  are  exerted  only  on  the  submerged  part  of  the 
airfoil  represented  by  masses  15  through  23  of  Figure  lb.  For  a  given  motion  input  p(t),  the  vector  FR  at  any 
time,  t,  at  or  after  the  step  is  given  by  : 


FRi=° 


,  i  =  1,  2..  .14 


pRi  (t)  =  JiL;  {p(t)  +  Lp(t)}  +  CNaLj 

2 


(  P(T)  +  p(x) )  r(t-T)dx 


,  i  =  15,  16... .23 


(5) 


where  CNa  is  the  static  normal  force  curve  slope  at  the  step  onset  angle,  Lj  is  the  span  of  the  ilh  airfoil 
element  (see  Fig.  lb),  and  T(t-x)  is  the  indicial  lift  response  given  by  either  Equation  (lb)  or  Equation  (lc). 
Notice  that  FR  is  not  a  function  of  the  generalized  coordinates  q^t),  q^t),  and  g(t)  since  in  the  ideal  response 
the  loading  is  given  by  (rigid  body)  2-D  airfoil  theory  alone.  Notice  also  that  if  p(t)  is  a  unit  step  function  of 
amplitude  A  a  and  L-  is  unity.  Equation  (5)  becomes  identical  to  Equation  (la)  (minus  the  initial  lift  C^Q). 

The  aeroelastic  force  vector  is  also  based  on  the  convolution  integral  formulation  for  an  airfoil  in  arbitrary 
motion,  and  again  acts  only  on  the  submerged  part  of  the  airfoil.  For  the  pitch  axis  at  the  quarter  chord  the 
loading  is: 


,  i  =  1,  2,... 14 
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F/U  (0  =  ^4  {-(<t»Nli  qi(t)  +  <j»N2i  CfcC1))  +  4>Ti  (  g(l)  +  ji(l) )}  + 


t  (6) 

CNaLi  I  [-(4>Kli  qi(t)  +  <t>N2i  q2<x))  +  <t>Ti  (gC*)  +  gOO)  ]  I>x)  dx 

Jo 

,  i  =  15,16,-23 

where  and  (j)^  are  the  values  of  the  NDOF1  and  NDOF2  mode  shapes  at  element  i,  respectively,  and 
is  the  value  of  the  TDOF  mode  shape.  The  negative  sign  on  the  NDOF  terms  is  due  to  the  sign  convention 
adopted.  The  non-integral  terms  are  the  so  called  apparent  mass  terms  since  they  may  be  moved  to  the  LHS  and 
combined  with  the  mass  term. 

The  inertial  loading  arises  from  the  fact  that  the  centroids  of  masses  13  through  23  do  not  coincide  with  the 
pitch  axis.  Thus,  angular  acceleration  about  the  pitch  axis  produces  an  inertial  normal  force  which  acts  at  the 
centroid.  The  inertial  force  vector  is: 


Fii=° 

Fjj  (t)  =  mj  r,  {  <(>Ti  g(t)  +  p(t)  > 


,i=  1.  2....12 
,  i  =  13,14,-23 


(7) 


where  nij  is  the  mass  of  element  i  and  rj  is  the  distance  from  the  pitch  axis  to  the  centroid. 

Aeroelastic  Moments 

The  moments  acting  on  the  inertial  elements  of  Figure  lb  are  also  decomposed  into  a  rigid  body  moment 
vector,  ,  an  aeroelastic  moment  vector,  TA  ,  and  an  inertial  moment  vector,  Tj  ,  giving  a  total  moment  of: 
T  =  Tj^  +  Ta  +  Tj .  The  rigid  body  moment  and  the  aeroelastic  moment  act  only  on  the  submerged  part  of  the 
airfoil.  The  expressions  for  these  components  are  simplified  by  the  fact  that  the  circulator}7  normal  force  acts  at 
the  quarter  chord  (for  a  flat  plate  airfoil)  giving  a  zero  moment  arm  in  the  present  study.  For  a  NACA  0015 
airfoil  the  circulator}7  normal  force  ( at  Re~  10^  )  acts  within  2%  fraction  of  chord  from  the  quarter  chord  and  on 
this  basis  has  been  neglected.  The  rigid  body  moment  for  rotation  about  the  quarter  chord  is: 


TRi  =  ° 

TRi  (t)  =  -JiLj  {  p(t)  +  —  p(t) } 
8 


,i  =  1.  2,-14 
,  i  =  15,16,-23 


(8) 


The  aeroelastic  moment  vector  acting  on  an  airfoil  in  arbitrary  motion  with  the  pitch  axis  at  the  quarter 
chord  is  given  by: 


TAi  =  ° 


,  i  =  1,  2,-14 


TAj  (t)  =  JiLj  {  qi(0  +  <t>N2i  q2(0)  -  <PTi  (  g(0  +|-g(t)) } 

,i  =  15,16,-23 

The  inertial  moment  vector  is: 

Tjj  =  JiP(t)  ,  i  =  1,2,... 12 


(9) 


12 


Tli  (t)  =  m;  r;  (<J>Nii  q^t)  +  <t>N2i  (fcW)  *  Ji  P(9 

,  i  =  13, 14,... .23 

where  Jj  is  the  mass  moment  of  inertia  of  element  i. 

Notice  that  Equations  (3a,b,c)  are  coupled  through  the  force  and  moment  terms  given  by  Equations  (5) 
through  (10). 

Input  Rigid  Body  Motion 

Two  input  motions  have  been  considered  and  are  illustrated  in  figure  4.  The  first  motion  is  a  step  function 
input  of  magnitude  A  a  arri  the  second  input  is  small  amplitude  ramp  modeled  after  the  experimental  angle  of 
attack  of  Figure  2.  These  motions  are  expressed  mathematically  by: 

p(t)  =  A  a  p(t-0)  (step)  (1  la) 

p(t)  =  4cl  t  (|i(t-0)  -  n(t-Ts))  +  A  a  p.( t-x  s)  (ramp)  (1  lb) 

Ts 

where  \i  is  the  unit  step  function,  and  xs  (=0. 13  semichords)  is  the  ramp  duration  w’hich  is  determined  from  the 
experimental  angle  of  attack  data.  As  will  be  seen,  it  will  be  necessary  to  take  the  Laplace  transforms  of 
Equations  (11a)  and  (1  lb)  which  yields,  respectively: 


P(s)  =  (step) 

$ 

(He) 

P(s)  =  A2L(1  -  e'STs)  (ramp) 

s2ts 

(lid) 

where  ,  s,  is  the  Laplace  variable. 

Solution  of  the  Generalized  Equations  of  Motion 

The  Laplace  transform  method  is  used  to  transform  the  three  generalized  differential  Equations  (3a,b,c)  into 
three  algebraic  equations  which  are  linear  in  terms  of  the  transformed  generalized  coordinates.  These  equations 
are  solved  simultaneously  for  the  transformed  generalized  coordinates  as  functions  of  the  Laplace  variable,  s. 
The  inverse  Laplace  transform  is  performed  to  transform  these  solutions  to  the  time  domain  to  obtain  q^(t), 
q^(t),  and  g(t).  The  Laplace  transform  of  each  of  Equations  (3)  can  be  done  by  hand  using  standard  tables.  To 
perform  the  extensive  algebra  required  in  the  solution,  the  computer  program  MACSYMA  has  been  used. 
MACSYMA  is  capable  of  symbolic  mathematics  required  in  the  solution  of  the  simultaneous  equations  in 
terms  of  the  Laplace  variable,  s.  The  resulting  equations  of  motion  in  the  Laplace  transform  domain  are  quite 
long  and  are  omitted  here.  Once  the  Laplace  transforms  of  Equations  (3)  have  been  done,  the  solution  proceeds 
by  first  rearranging,  using  MACSYMA,  the  system  of  three  equations  into  the  following  form: 
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An(s)  A12(s)  A13(s) 
A 2  i(s)  A22(s)  A23(s) 
A3i(s)  A32(s)  A33(s) 


rQi(s)' 

fFl(s)\ 

1 

ll 

Q2(S) 

l  =  i 

f  ) 

F2(s) 

1 

l  G(s)  J 

i  1 

\F3(s)  ) 

(12) 


whereQi(s)  is  the  Laplace  transform  of  q-j(t),  Q2(s)  of  q2(t),  and  G(s)  of  g(t).  In  Equation  (12),  P(s)  is  the 
Laplace  transform  of  the  motion  input  parameter  p(t)  (Equation  (11c)  or  (1  Id)).  Equation  (12)  can  be  solved 
using  Cramer's  rule  to  yield: 


(X(s)}  =  (R(s)}  P(s) 


(13a) 


where  {X(s)>  is  a  3x1  vector  given  by: 

|Ql(s)| 
w4(s) 
\G(s)  ] 


(13b) 


and  the  3x1  vector  (R(s)}  is  an  aeroelastic  transfer  function  which  is  a  function  of  the  Laplace  variable  ,  s,  only 
and  is  given  by: 


Rl(s)= 


Fi(s) 

AU(s) 

Ais(s) 

Ai  i(s) 

Fi(s) 

Ai3(s) 

Ai  i(s) 

Ai2(s) 

Fi(s) 

F2(s) 

A22(s) 

A23(s) 

A2  i(s) 

F2(s) 

A23(s) 

A2  i(s) 

A22(s) 

F2(s) 

F3(s) 

A32(s) 

A33(s) 

A3i(s) 

F3(s) 

A33(s) 

-  ,  R^(s)= 

A3  i(s) 

A3  2(s) 

F3(s) 

An(s) 

Ai  2(s) 

Ai3(s) 

An(s) 

Ai2(s) 

Ai3(s) 

Ai  i(s) 

Ai  2(s) 

Ai3(s) 

A2 1  (s) 

A22(s) 

A23(s) 

A2  i(s) 

A22(s) 

A23(s) 

A2  i(s) 

A22(s) 

A23(s) 

A3  i(s) 

A32(s) 

A33(s) 

A3  i(s) 

A3  2(s) 

A33(s) 

A3  l(s) 

A3  2(s) 

A33(s) 

where  the  brackets  indicate  the  determinant  It  has  been  found  efficacious  to  perform  the  necessary  algebra  so 
that  Rj(s),  R2(s),  and  R3(s)  are  the  form  of  the  quotients  of  two  polynomials.  This  amounts  to  manipulating 
Equation  (12)  so  that  all  of  the  terms  in  the  square  matrix  [A]  and  the  terms  in  the  vector  {F}  are  at  most 
polynomials  (i.e.  not  themselves  quotients).  Quotients  will  appear  in  [A]  and  {F}  due  to  the  exponential  form 
of  theindicial  function  T(t-t).  Recalling  form  Equations  (lb)  and  (lc),  the  indicial  function  is  of  the  form. 

F(t)  =  1  -  Aeat  -  Be^1 ,  which  has  the  Laplace  transform: 

M-1  A.  B  -  (l-A-B)s2  -f  (Ab-t-Ba-a-b)s  +ab 
^  s  s-a  s_b  s(s-a)(s-b) 


(13c) 


The  quotient  may  be  eliminated  by  simply  multiplying  the  entire  system  of  Equations  (12)  by  the  denominator 
of  y(s).  This  will  render  the  terms  of  [A]  and  {F}  as  non-quotients  and  simplify  the  Laplace  inversion  by 
partial  fraction  decomposition  which  follows. 

Substitution  of  a  particular  input  motion,  P(s),  into  Equation  (13a)  gives  the  solution  for  the  generalized 
coordinates  in  the  Laplace  domain.  For  the  step  input  of  Equation  (11c),  the  result  for  each  generalized 
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coordinate  is  in  the  form  of  a  quotient  of  two  polynomials  in  s,  wherein  the  numerator  is  a  tenth  order 
polynomial  and  the  denominator  is  an  eleventh  order  polynomial.  For  the  ramp  motion  of  Equation  (1  Id)  the 
result  is  a  tenth  order  divided  by  a  twelfth  order  polynomial. 

In  either  case,  the  inverse  Laplace  transform  is  accomplished  by  decomposing  the  quotients  using  partial 
fractions  into  functions  which  can  be  inverted  by  hand.  In  the  partial  fraction  decomposition  MACSYMA  has 
been  used  to  find  the  roots  of  the  polynomial  in  the  denominator.  In  the  step  input  case  the  denominator  has 
five  real  roots  while  in  the  ramp  input  the  denominator  has  six  real  roots  with  a  repeated  root  at  s=0.  The  real 
roots  arise  from  the  leading  constant  and  the  exponential  terms  in  the  indicial  response  function,  T,  and  the  real 
roots  of  the  motion  function  P(s).  In  both  the  step  and  ramp  input  cases,  the  denominator  further  has  three 
pairs  of  complex  conjugate  roots  which  represent  the  oscillatory  (harmonic)  response  of  the  structure  at  the  three 
frequencies  associated  with  the  NDOF1,  NDOF2,  and  TDOF.  The  numerical  values  of  the  frequencies  of  the 
harmonic  response  can  be  traced  to  the  mass  (including  apparent  mass)  and  stiffness  terms  in  the  governing 
aeroelastic  equations.  The  set  of  simultaneous  equations  for  the  undetermined  constants  in  the  partial  fraction 
expansion  terms  have  been  solved  using  direct  factorization  with  maximal  column  pivoting  [10].  Pivoting  is 
recommended  since  the  coefficients  in  the  simultaneous  equations  vary  over  several  orders  of  magnitude. 

It  is  worth  noting  that  in  the  ramp  motion  input  case  it  is  not  necessary  to  use  the  full  form  of  P(s)  as 
given  by  Equation  (lid).  This  is  because  the  last  exponential  term  on  the  right  hand  side  results  in  an  inverse 
transform  identical  to  the  inverse  transform  of  the  term  which  precedes  it,  but  shifted  in  time  by  an  amount,  x§. 
Thus  it  is  only  necessary  to  use  P(s)  =  Aa/s^x§  when  calculating  this  case,  and  subsequently  subtracting  from 
this  solution  an  identical  solution  shifted  in  time  by  xs-  This  effect  can  be  seen  in  the  solution  presented  below 
as  Equation  (14b). 

For  the  step  input  case,  the  solutions  for  the  generalized  coordinates  q^(t),  and  g(t)  can  be  represented 
generally  by  the  function  f(t)  which  has  the  following  form: 

f(t)  =  A  +  Bebt  +  Cect  +  Dedt  +  Feft  + 

eht(Gcoso>it  +  Hsinwit)  +  e)  t(Jcoso)2t  +  Ksino^O  +  emt(Mcoso)3t  +  Nsino^t) 

where  the  constants  are  given  in  the  following  table.  The  constants  in  Tables  2  and  3  were  computed  using  the 
Wagner  function  for  the  indicial  lift  function,  T. 

Table  2.  Constants  for  Generalized  Coordinates  for  Step  Input  using  Wagner  Function 


qx(t) 

<fe(0 

g(0 

A 

1.241e-2 

1.196e-4 

0 

B 

-1.417e-3 

-3.178e-6 

-9.786e-9 

C 

-6.001e4 

-1.537e-5 

7.451e-9 

D 

-3.162e-3 

1.601e-5 

-1.531e-7 

F 

-3.619e-5 

4.268e-5 

-2.328e-10 

G 

8.407e-3 

-3.149e-5 

-5.536-5 

H 

2.042e-2 

-9.484e-5 

4.836e-5 
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J 

7.314e-4 

7.155e-3 

-3.273e4 

K 

1.476e-4 

1.725e-3 

-5.791e-5 

M 

-1.584e-2 

-6.993e-3 

-2.075e-2 

N 

-3.740e4 

-1.377e-4 

6.344e4 

The  remaining  constants  apply  to  each  of  qj(t),  (£>(t),  and  g(t) 

b  =  -0.04522 

c  =  -0.04565 

d=  -0.2877 

f= -0.29989 

h  =  -0.55591 

j  =  -0.49752 

m  =  -2.71184 

co  1(NDOFl)  =  2.3361 

o)2(NDOF2)  =  12.186 

co3(TDOF)  =  52.2147 


For  the  ramp  motion  input,  the  solutions  for  q  j(t),  q2(t),  and  g(t)  in  the  time  domain  can  be  represented  by 
a  function  of  the  form:  f(t)  -  p(t-xs)  f(t-xg)  (e.g.  q^t)  =  f(t)  -  p(t-xs)  f(t-xs)  ),  where  xs  is  the  ramp  duration 
(0.13  semichords),  and  p(t-xs)  is  the  unit  step  function.  As  has  been  mentioned,  the  shifting  in  time  of  the 
function  f(t-xs)  is  due  to  the  exponential  term  in  the  Laplace  transform  of  the  motion  input  P(s)  of  Equation 
(lid).  In  the  ramp  case,  f(t)  has  the  form: 

f(t)  =  A  +Et  +  Bebt  +  Cect  +  Ded'  +  Feft  +  (J4b) 

e^1  '(Gcosco  j  t  +  Hsincoit)  +  d '(loosest  +  Ksinco2t)  +  e™  *(Mcosa>  3 1  +  Nsinco3t) 

Notice  the  ramp  motion  gives  rise  to  a  term  proportional  to  time  t.  This  might  be  expected  since  the  ramp 
motion  produces  an  upwash  at  the  3/4  chord  proportional  to  t.  For  times  beyond  xg,  this  term  goes  to  Exs 
(=const)  to  give  the  steady  state  lift  The  constants  in  f(t)  and  are  given  in  the  following  table. 


Table  3.  Constants  for  Generalized  Coordinates  for  Sudden  Ramp  Input  using  Wagner  Function 


qi(0 

g(0 

A 

-0.3578 

-2.91e-3 

-2.874e4 

B 

0.3246 

-1.335e-3 

8.305e-7 

C 

1.759e-2 

4.467e-3 

4.240e-7 

D 

8.482e-2 

4.460e4 

3.971e-6 

E 

9.548e-2 

9.205e4 

0 

F 

5.894e4 

1.114e-3 

1.20e-7 

G 

-6.988e-2 

3.190e-5 

1.918e4 

H 

1.106e-2 

-2.777e-5 

-1.367e4 
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J 

-1.135e4 

-1.278e-3 

4.527e-5 

K 

4.573e-4 

4.463e-3 

-2.046e4 

M 

1.644e-4 

6.946e-5 

4.579e-5 

N 

-2.327e-3 

-1.027e-3 

-3.055e-3 

*The  values  of  the  remaining  constants  are  the  same  as  given  in  Table  2. 

The  above  solutions  are  in  closed  form  and  have  the  mathematical  properties  that  might  be  anticipated. 
That  is,  the  exponential  terms  giving  the  response  to  a  circulatory  build-up  in  lift  after  the  motion  input,  and  an 
oscillatory  component  which  occurs  at  three  distinct  frequencies  that  arise  due  to  the  modal  excitation  and 
coupling  in  the  forcing  terms  on  the  RHS  of  Equations  (3).  Notice  the  reduction  in  the  aeroelastic  frequencies 
in  the  results  above  from  the  natural  frequencies  (given  previously)  in  the  respective  modes  due  to  apparent  mass 
effects.  The  solutions  given  above  are  plotted  in  Figures  5a-l  where  the  modal  coupling  can  be  observed  in  the 
frequency  content  of  each  modal  solution.  The  solutions  are  given  for  both  the  step  and  ramp  motions  and 
have  been  calculated  using  both  the  Wagner  function  and  the  Kussner  function  as  the  indicial  lift.  The 
conditions  are  indicated  in  the  caption  of  each  figure.  As  might  be  expected  the  step  input  tends  to  excite  the 
higher  frequency  TDOF  more  than  the  more  gradual  ramp  input. 

A  comparison  of  Figures  5g  and  5h  as  well  as  Figures  5i  and  5j  indicates  that  the  oscillations  in  the 
Kussner  function  solution  persist  for  a  longer  period  of  time  than  the  Wagner  solution.  This  is  due  to  the  fact 
that  the  Kussner  function  provides  no  damping  at  the  step  onset  where  the  Kussner  function  has  a  value  of  zero. 
The  Wagner  function,  on  the  other  hand,  has  an  instantaneous  value  of  0.5  at  the  step  origin  which  provides 
significant  damping  through  the  integral  terms  in  Equations  (6). 

It  is  of  course  possible  to  input  other  motions  than  those  of  Equations  (11a)  and  (1  lb).  For  example,  the 
trigonometric  function:  p(t)  =  (Aa^jiX^s1  "  sincostX[i(t-0)  -  |x(t-xs))  +  Aa[x(t-xs)  may  be  used  to  represent 
the  step,  where  cos  =  2kIx$.  This  function  has  continuous  first  and  second  derivatives.  The  solution  for  the 
generalized  coordinate  is  of  the  form  f(t)  -  |fit-xs)  f(t-xg) ,  where: 

fit)  =  A  +Et+(Usinoost  +  Vcoscost)  +  Be'bt  +  Ce’*  +  De‘dt  +  Fe^  + 
e~ht(Gsino) \t  +  Hcosoo  it)  +  e‘Jt(Jsinco2t  +  Kcosc02t)  +  e~mt(Msino)3t  +  Ncoso)3t) 

Notice  that  the  period  of  cos  is  xg  and  consequently,  the  harmonic  terms  containing  cog  cancel  identically  for 
times  greater  thanxg  ( according  to  fit)  -  ^(t-xg)  f(t-xg)  ).  The  results  for  this  input  are  similar  to  the  ramp 
input  and  are  not  presented  here. 

Calculation  of  the  Sensible  Normal  Force  at  the  Load  Cell 

The  purpose  of  the  present  analysis  is  to  determine  the  effects  of  aeroelasticity  on  the  output  of  the  strain 
gauge  load  cell  used  in  the  experiments  described  above.  In  these  experiments  the  strain  gauge  output  is 
interpreted  ideally  as  that  due  to  a  moment  exerted  by  a  normal  force  applied  at  the  midspan  of  the  submerged 
part  of  the  airfoil  (point  loads  were  applied  here  to  calibrate  the  strain  gauge  bridge).  In  reality,  however,  this 
is  not  the  case  since  the  strain  gauge  output  is  determined  solely  by  the  instantaneous  beam  curvature  at  the  cell. 
This  curvature  is  due,  not  to  aerodynamic  loading  alone,  but  rather  the  total  aeroelastic  structural  response. 
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Figures  5c  and  5d.  Step  Solution  for  the  NDOF2  Generalized  Coordinate  q2(t)  using: 
a)  Wagner  Function,  b)Kussner  Function. 
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Figures  5e  and  5f.  Step  Solution  for  the  TDOF  Generalized  Coordinate  g(t)  using: 
a)Wagner  Function,  b)Kussner  Function. 
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Figures  5g  and  5h.  Ramp  Solution  for  the  NDOF1  Generalized  Coordinate  qj(t)  using: 
a)  Wagner  Function,  b)Kussner  Function. 
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Figures  5i  and  5j.  Ramp  Solution  for  the  NDOF2  Generalized  Coordinate  q2(t)  using: 
a)  Wagner  Function,  b)Kussner  Function. 
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Figures  5k  and  51.  Ramp  Solution  for  the  TDOF  Generalized  Coordinate  g(t)  using: 
a)  Wagner  Function,  b)Kussner  Function. 
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From  beam  theory  the  moment  at  the  center  of  the  load  cell  (centroid  of  element  7)  is  related  to  the  beam 
curvature  by: 

MLdt)  =  EI{qi(t)^L>)  +q2(t)^^.)  >  (15a) 

dx2  ' 7  dx2 

where  x  is  measured  along  the  span.  The  second  derivative  has  been  computed  numerically  from  the  NDOF 
mode  shapes  for  elements  5  through  9  using  a  five  point  numerical  derivative  given  by  Richardson's 
extrapolation  method  [10].  A  time  dependent  "sensible”  normal  force  acting  at  the  midspan  of  the  submerged 
portion  of  the  airfoil  which  produces  the  same  moment  as  that  given  by  Equation  (15a)  may  then  be  defined  for 
comparison  with  experimental  normal  force  data.  For  this  purpose,  the  moment  arm  between  the  mass  7  and 
the  midspan  is  L =  37.0  in.  These  results  are  presented  in  the  form  of  a  sensible  normal  force  coefficient 
definedas: 


Qsift)  =  CNo  +  MLC(t)/LLC  (15b) 

pU2bL 

w'here  is  the  normal  force  at  the  origin  of  the  motion,  p  is  the  density  of  w7ater,  and  L  is  the  submerged 
airfoil  length  (42.0  in),  U  and  b  are  the  freestream  velocity  and  the  semichord  length,  respectively. 

Aeroelastic  Analysis  Results 

Sensible  force  calculations  using  the  present  model  are  compared  below  with  the  experimental  normal  force 
data  of  Reference  [4].  A  comparison  is  also  made  with  recent  accelerometer  data  taken  on  the  Ohio  U.  tow7  tank 
rig  undergoing  the  motion  illustrated  in  Figure  2.  In  these  tests  an  accelerometer  (PCB  Flexcel  Series  336A) 
w?as  placed  on  mass  13  (see  Figure  lb)  with  the  accelerometer  centered  on  the  pitch  axis.  The  accelerometer 
data  have  been  integrated  once  and  put  in  the  form  of  velocity.  The  present  model  can  predict  velocity  at  any 
point  on  the  structure  by  taking  the  time  derivative  of  Equations  (2).  The  accelerometer  data  and  normal  force 
data  below  were  acquired  in  separate  runs. 

Sensible  Force  and  Accelerometer  Comparisons 

The  two  motion  inputs  given  by  Equations  (11a)  and  (11b)  have  been  considered.  The  sensible  force 
corresponding  to  the  step  input  is  compared  with  experimental  strain  gauge  data  in  Figures  6a  and  6b.  Figure  6a 
has  been  computed  using  Wagner's  function  and  Figure  6b  has  been  computed  using  Kussner's  function.  .Also 
shown  is  the  theoretical  response  based  on  the  Wagner  function  and  Kussner  function,  respectively.  There  are 
three  frequencies  present  in  the  aeroelastic  model.  The  lowest  frequency  has  a  period  of  approximately  2.5 
semi  chords  and  is  associated  with  the  NDOF1  fundamental  mode.  The  second  frequency  has  a  period  of  about 
0.5  semichords  and  is  due  to  the  NDOF2  mode,  and  the  third  highest  frequency  with  a  period  near  0. 1  semichord 
is  due  to  the  coupling  between  the  TDOF  and  the  NDOF1  and  NDOF2.  It  is  clear  that  aeroelastic  reactions 
have  caused  significant  deviation  from  the  theoretical  response  particularly  before  about  1.5  semichords. 

The  sensible  force  based  on  the  aeroelastic  response  with  the  sudden  ramp  input  of  Equation  (lib)  is 
illustrated  in  Figures  6c  and  6d.  The  response  directly  after  the  motion  inception  is  in  better  agreement  with 


u.  /u 


a)  Wagner 


Figures  6a  and  6b.  Sensible  Force  Results  for  Step  Input  using 
a)  Wagner  function,  b)  Kussner  function 


Figures  6c  and  6d.  Sensible  Force  Results  for  Ramp  Input  using 
a)  Wagner  function,  b)  Kussner  function 
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the  strain  gauge  data  than  the  step  input  results.  A  comparison  of  the  aeroelastic  response  of  Figures  6a,b  and 
6c,d  indicates  that  the  step  input  excites  the  NDOF  and  TDOF  coupling  more  than  the  sudden  ramp  motion. 
At  small  times  the  experimental  data  exhibit  discernible  oscillations,  which  in  the  analytical  force,  are  associated 
with  the  NDOF2  mode.  The  amplitude  of  these  oscillations  is  in  better  agreement  with  the  Wagner  function 
model  than  the  Kussner.  As  has  been  pointed  out  in  relation  to  Figures  5i  and  5j,  the  Kussner  function 
provides  very  little  damping  of  the  NDOF2  mode,  which  is  also  evident  in  the  sensible  force  results  of  Figure 
6d.  It  appears  that  the  use  of  the  Wagner  function  gives  better  agreement  overall  with  the  experimental  data. 
It  should  be  noted  that  the  uncertainty  in  the  experimental  data  is  near  10%. 

A  comparison  between  the  aeroelastic  analysis  using  motion  1  step  input  and  accelerometer  data  is  shown 
in  Figure  7a  (based  on  Wagner  function).  The  accelerometer  data  have  been  integrated  numerically  and  put  in 
the  form  of  velocity  data.  Again,  the  use  of  the  step  function  input  excites  the  TDOF  mode  more  than  the 
experimental  data  indicates.  The  use  of  the  sudden  ramp  motion  gives  the  results  of  Figure  7b  (based  on 
Wagner  function)  where  the  agreement  between  the  experimental  data  and  the  analysis  is  very  good.  Use  of 
Kussner's  function  gives  results  similar  to  those  in  Figures  7a  and  7b. 

Part  II:  Inverse  Aeroelastic  Analysis 

Described  below  is  an  inverse  aeroelastic  analysis  directed  towards  computing  the  aerodynamic  indicial 
response  given  the  system  structural  properties  and  the  total  aeroelastic  response.  The  aeroelastic  response 
may  be  known  from  experimental  strain  gauge  and/or  accelerometer  data  which  by  their  very  nature  will 
contain  the  total  aeroelastic  response.  The  following  analysis  focuses  on  computing  Wagner’s  function 
from  the  analytical  aeroelastic  response  which  has  been  described  in  detail  above.  Only  the  ramp  input 
motion  is  considered. 

The  problem  in  the  inverse  analysis  may  be  stated  as  follows:  Given  the  structural  properties  of  a 
system(i.e.  mass  and  stiffness)  and  the  time  dependent  aeroelastic  response  to  a  known  excitation  (motion 
input),  what  must  the  aerodynamic  indicial  response  be  to  produce  this  aeroelastic  response? 


Indicial  Response  Formulation 

The  Laplace  transform  method  has  been  used  in  the  present  inverse  analysis.  Taking  the  Laplace 
transform  of  Equation  (3a),  simplifying,  and  rearranging  yields  an  expression  of  the  form: 


(ais  +  — )Qi(s)  +  a3sQ2(s)  +  (214  +a5s)(P(s)  +  G(s))  = 
s 

[a6sQi(s)  +  a7sQ2(s)  +  ag(l  +  s)(P(s)  +  G(s))]y  (s) 


(16) 


where  ,  s,  is  the  Laplace  variable,  a1  through  a^  are  constants,  and  as  before,  Qj(s),  Q2(s),  and  G(s)  are  the 
Laplace  transforms  of  the  time  dependent  generalized  coordinates  q^(t),  q2©,  and  g(t).  The  term  P(s)  is  the 
Laplace  transform  of  the  ramp  motion  input  (Equation  (lid)  given  by:  P(s)  =  (1  -  e"sxs)  ,  y  (s)  is  the 

Laplace  transform  of  the  indicial  response  function  (Wagner  function)  which  is  to  be  determined.  It 
should  be  noted  that  the  terms  in  Equation  (16)  containing  P(s)  and  G(s)  have  been  grouped  together.  This 
is  a  due  to  the  particular  structure  being  modeled  wherein  the  value  of  the  mode  shape  in  the  TDOF  is  unity 
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Figures  7a  and  7b.  Comparison  of  Velocity  Predictions  (using  Wagner  function)  and 
Integrated  Accelerometer  Data  for  Step  and  Ramp  Input  Motions 
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for  elements  13  through  23  (see  Figures  lb  and  3).  The  result  is  that  the  terms  containing  the  RDOF 
coordinate  P(s)  and  the  TDOF  coordinate  G(s)  have  the  same  coefficients. 

Equation  (16)  may  further  be  simplified  by  eliminating  the  TDOF  variable  G(s)  as  follows.  Taking 
the  Laplace  transform  of  the  TDOF  Equation  (3c)  and  solving  for  G(s)  yields  an  expression  of  the  form: 

G(s)  =  ( - 1 - )((c4s  +  C5S2)P(s)  +  cgs2Qi(s)  +  C7S2Q2(s))  (17) 

c^s2  +  C2S  +  C3 

where  c^  through  c^  are  constants.  Notice  that  G(s)  is  not  a  function  of  the  indicial  response  y  (s)  which 
is  a  consequence  of  the  pitch  axis  being  located  at  the  quarter  chord.  From  an  experimental  design 
standpoint  this  is  a  judicious  choice  since  it  simplifies  the  equations  of  aeroelasticity.  Substituting 
Equation  (17)  into  (16),  simplifying  to  eliminate  any  superfluous  constants,  and  solving  for  the  indicial 
function  y  (s)  yields: 

(.)  -  Bi(s)Qi(s)  +  B2(s)Q2(s)  +  B3(s)P(s)  (18) 

B4(s)Qi  (s)  +  B5(s)Q2(s)  +  B6(s)P(s) 

where: 

Bi(s)  =  bi  +SL+-  sfe+.b4.s) .  (19a) 

s2  s2  +  b5S  +  bg 

B2(s)  =  l>7  +  -.s^tM..  (19b) 

s2  +  b5S  +  bg 

Bg  (s)  +  bi  i  +  (bl2  +  bi3s+  bi4s2)  (19c) 

s  s2  +  b5S  +  b6 

B4(s)  =  bi  5  +  ■  bl  6S(1+S)  (19d) 

s2  +  b5S  +  b6 

B5(s)  =  bi  7  +  blgs(1+s)  (19e) 

s2  +  b5S  +  b$ 

Be(s)  =  b!  9(1+1)  +  <b20  +  b2is+b22S2.)  (190 

s  s2  +  b5S  +  b6 

where  the  lower  case  b1  through  b22  are  known  constants,  given  in  Appendix  A.  The  physical  origins  of 
the  constants  b ^  through  b22  may  of  course  be  traced  to  the  constants  in  the  Equations  (3a)  and  (3c).  For 
example,  b2  is  the  generalized  stiffness  Kj  in  Equation  (3a).  Variations  of  Equations  (19)  are  of  course 
possible  depending  on  how  one  chooses  to  perform  the  simplifying  algebra.  The  forms  above  have  been 
found  efficacious  in  that,  as  will  be  seen,  it  will  be  necessary  to  express  the  Laplace  variable  ,  s,  as  a 
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complex  number  and  then  separate  the  resulting  real  and  imaginary  parts  of  the  B^(s)  through  Bg(s) 
functions  of  Equations  (19). 

The  principal  task  is  to  perform  the  inverse  Laplace  transform  of  Equation  (18)  for  known  Q^(s), 
Q2(s),  and  P(s).  In  situations  where  qj(t),  q2©  ,  and  p(t)  are  known  analytical  functions  of  time  it  will 
generally  be  possible  to  express  Qj(s),  Q2(s),  ^d  P(s)  as  functions  of  s  by  simply  taking  their  Laplace 
transforms.  In  this  case  it  would  theoretically  be  possible  to  substitute  the  transformed  functions  into 
Equation  (18)  and  invert  the  resulting  expression  for  y  (s)  analytically  using  partial  fraction  decomposition. 
This  is  essentially  the  approach  that  has  been  used  in  the  aeroelastic  analysis  described  in  the  Part  I  of  this 
report,  where  the  transformed  indicial  function  y  (s)  has  been  expressed  in  analytical  form  and  the  resulting 
expressions  for  Qi(s),  Q2(s),  and  G(s)  have  been  inverted  analytically  for  a  given  motion  input  P(s). 

It  is  a  purpose  of  this  research  to  investigate  a  method  for  experimentally  determining  the  indicial 
function  and,  as  such,  it  is  expected  that  the  functions  qj(t),  q2(t),  the  motion  input  p(t)  will  be  in  the 
form  of  discrete  data  in  time  which  generally  cannot  be  curve  fit  with  simple  closed  form  functions.  As 
will  be  seen,  in  this  instance  it  is  possible  to  perform  a  numerical  Laplace  inversion  of  Equation  (18).  As 
discussed  below,  the  functions  q^(t)  and  q2(0  may  be  obtained  from  the  total  aeroelastic  response  contained 
in  strain  gauge  and/or  accelerometer  data. 

In  the  above  analysis  Equation  (3a)  (NDOF1)  and  (3c)  (TDOF)  have  been  combined  to  give  Equation 
(18).  It  is  alternatively  possible  to  use  Equations  (3b)  (NDOF2)  and  (3c),  however,  it  has  been  found  that 
the  use  of  (3a)  rather  than  (3b)  is  a  more  robust  formulation.  Although  no  systematic  error  analysis  has 
been  performed,  it  appears  that  determining  (from  experimental  data)  the  NDOF2  coordinate,  q2(t),  is  more 
susceptible  to  error  than  determination  of  the  larger  NDOF1  coordinate,  qj(t).  Equation  (3a)  is  less 
sensitive  to  these  errors  than  Equation  (3b)  since  the  coefficients  of  q2(0  (3a)  are  an  order  of  magnitude 

smaller  than  those  in  (3b).  It  is  also  possible  to  combine  all  three  Equations  (3a),  (3b),  and  (3c).  This 
method  has  likewise  been  found  to  be  less  stable  than  that  based  on  Equations  (3a)  and  (3c). 

Numerical  Inverse  of  Laplace  Transform 

The  method  of  Reference[l  1]  has  been  used  to  numerically  invert  Equation  (18).  Briefly,  the  method 
evaluates  the  inverse  by  means  of  a  Fourier  series,  and  is  essentially  a  numerical  computation  of  the 
Laplace  inversion  formula: 

i-  a+ioo 

T(t)  =  -J—  esty(s)ds  (20) 

2jli  I  -nr* 

J  a-ioo 

where  ,  a,  is  a  constant  such  that  any  singularities  of  y  (s)  lie  to  the  left  of  the  line  integral  in  the  complex 
plane.  Letting  s  =  a  +  ico  so  that  y  (s)  =  y(a,co) ,  Equation  (20)  may  alternatively  be  expressed  by: 


T(t)  =  -  f  [  Re  ( y  (a,co) )  cosrnt  -  Im  ( y  (a,co ) )  sincot  ]  dco 

Jt  i  0 


(21) 


A  trapezoidal  rule  approximation  to  Equation  (21)  on  the  time  interval  (0,2T)  is  given  by: 
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T(t)  =  ( ^1)  f  iM  +  j  {  Re  ( y  (a ,  fcrc/T) )  cos(toit/T)  -  Im  ( y  (a  ,  kjt/T) )  sin(k:Jit/T)  J  ]  (22) 

T  2  k=l 

Equation  (22)  provides  an  approximation  to  the  exact  value  of  the  indicial  response  function.  The  value  of 
T  is  chosen  to  be  T  =  0.8  t,^  ,  where  is  the  largest  time  for  which  T(t)  is  to  be  approximated. 
The  value  of  a  is  chosen  according  to  a  =  a  -  ln(E)/2T  ,  where  a  is  a  number  slightly  larger  than  the 
maximum  real  pole  of  y  (s),  and  E  is  the  desired  maximum  relative  error.  The  relative  error  is  defined  as 
the  error  in  the  approximation  divided  by  the  maximum  value  of  the  function  to  be  approximated.  Thus,  if 
the  approximation  is  to  converge  to  at  least  two  significant  figures  the  relative  error  would  be  set  to  0.005. 
Of  course  in  practice  the  summation  in  Equation  (22)  would  contain  only  a  finite  number  of  terms 
designated  henceforth  as  N.  The  reader  may  consult  Reference  [1 1]  for  additional  details  of  the  inversion 
method. 

To  numerically  invert  Equation  (18)  using  Equation  (22),  it  is  necessary  to  treat  the  Laplace  variable,  s, 
as  a  complex  number  and  separate  y  (s)  into  real  and  imaginary  parts.  The  functions  B  |  (s)  through  Bg(s) 
as  well  as  the  transformed  generalized  coordinates  Qj(s)  and  C?2(s)  and  the  motion  input  P(s)  will  each 
contain  a  real  and  imaginary  part.  Letting  s  =  x  +  iy,  the  resulting  expression  for  y  (s)  will  be  given  by: 

Y(x,y)  =  YR(x,y)  +  iyi(x,y)  t233) 

where  y  R(x,y)  and  y  j(x,y)  are  the  real  and  imaginary'  parts,  respectively,  of  y(x,y)  which  are  given  by: 

YR(*,y)  =  NRSR*NlSr.  ,  Yl(x,y)  =  NlSR'NHS-l  (23b) 

Sr+S?  sR+s? 

and: 

Nr(x,)0  =  (  B1RQ1R  -  BnQn  )  +  (  B2rQ2R  -  B21Q21  )  +  (  B3RPR  -  )  (23c) 

Ni(x,y)  =  (  BirQh  +  BhQir  )  +  (  B2rQ2I  +  B2iQ2r  )  +  (  B3RP1  +  B31PR )  (23d) 

SR(x,y)  =  ( B4RQ1R  -  B41Q11 )  +  (  B5RQ2R  -  B51Q21 )  +  (  B6rPr  -  B61P1 )  (23e) 

Si(x,y)  =  (  B4rQh  +  B4IQ1R  )  +  (  B5RQ21  +  B51Q2R  )  +  (  B6rPi  +  B6IPr  )  (23f) 

where  the  subscripts  R  and  I  indicate  real  and  imaginary  parts,  respectively.  The  separation  of  the 
generalized  coordinates  and  the  motion  input  into  real  and  imaginary  parts  is  discussed  below.  The 
separation  of  the  functions  Bj(s)  through  B^(s)  is  accomplished  by  substituting  s  =  x  +  iy  into  the 
Equations  (19)  and  performing  the  necessary  algebra.  For  example,  the  function  B3(x,y)  =  B3p(x,y)  + 
iB3i(x,y),  has  real  and  imaginary  parts  given  by: 


31 


B3R(x,y)  =  bl°x  +  bi  i  +  (*2-y2+b5*+  b6Xbi2+bi3X+bi4(x2-y2))  +  (2xy+b5y)(bi3y+2bi4xy) 
x2+y2  (x2-y2+b5X+bg)2  +  (2xy+b5y)2 


(24) 


B3j(x  y)  —  ~  bl Oy  I  (x2-y2+b5X+b6Xbi3y+2bi4xy)  -  (2xy+b5y)(b12+bi3x+b14(x2-y2)) 
x2+y2  (x2-y2+b5X+b6)2  +  (2xy+b5y)2 

The  remaining  expressions  for  B  j(s)  through  B^(s)  may  similarly  be  separated  into  real  and  imaginary  parts. 
For  using  Equation  (24)  in  the  approximating  Equation  (22)  one  would  set  x  =  a  ,  and  y  =  kjtt/T.  It  should 
be  noted  that  at  points  of  discontinuity,  the  above  method  will  converge  to  T(t)  =  (  T(t+)  +  F(t-)  )/ 2  , 
where  t  is  the  time  of  the  discontinuity.  Indicial  responses  may  contain  points  of  discontinuity  such  as 
Wagner's  function  which  is  discontinuous  at  t  =  0. 

Other  numerical  inversion  methods  may  also  be  used  such  as  that  described  in  Reference  [12],  in  which 
the  inverse  is  obtained  in  terms  of  orthonormal  Laguerre  polynomials.  The  inverse  transform  is  given 
by: 

r(t)  =  «?*  2  Ch*n(i)  (2s) 

n=0  T 

where  ,  a,  is  a  real  constant  larger  than  any  real  singularities  of  y  (s),  and  <£n(l.)  are  Laguerre  polynomials. 

The  constant  T  is  a  scaling  factor  chosen  according  to:  T=  tmax/N  ,  w7here  ^  defines  the  time  interval  of 
interest,  and  N  is  the  number  of  terms  in  the  summation  of  Equation  (25).  This  choice  of  T  ensures  that 
the  oscillatory  frequencies  of  the  Laguerre  polynomials  are  compatible  w7ith  those  of  the  function  to  be 
approximated  The  constants  cQ  are  given  by: 

Cq  =  — L-j  h(0j)  ,  <^=—2—  2  h(0j)  cos  (n0j)  ,  n*0  (26a) 

N+l  j=0  N+l  j_Q 


where:  0;  =  (BL±_L)  2L  j  =  1, 2.....N  (26b) 

J  2N+1  2 

Q.  0.  0. 

h(0i)  =  J—  [  y  p(  a,  -1-cot—)  -  (J-cot-i.)  y  (  a,  J-cotJ-)  ]  (26c) 

J  2T  2T  2  2T  2  j  2T  2 

where  y  R(x,y)  and  y  j(x,y)  are  defined  above.  In  evaluating  Equation  (25)  using  Equations  (26a)-(26c)  one 

0. 

would  set  x  =  a,  and  y  =  J— cot— .  The  Laguerre  polynomials  in  Equation  (25)  may  be  evaluated  using  the 
2T  2 

recursion  formulas: 

$o(t)  =  e-t/2 

<b1(t)  =  (l  -  t)<L>0(t) 


(26d) 
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n®n(t)  =  (2n-l-t)*n.1(t)-(n-l)*n-2(t)  ,  n>l 

It  has  been  found  that  this  method  agrees  with  the  method  of  Reference  [12]  for  intermediate  and  large 
values  of  time,  however,  for  small  time  the  method  gives  oscillatory  results.  Results  from  both  methods 
are  presented  below. 

Transform  of  Generalized  Coordinates:  Analytical  Form 

The  numerical  inversion  of  Equation  (18)  has  been  done  using  the  analytical  solutions  for  qj(t)  ,q2(t)  » 
and  g(t)  given  by  Equation  (14b)  for  the  ramp  motion  input  case.  In  this  case  the  functions  Qj(s)  and 
Q2(s)  are  in  the  form:  F(s)  (1  -  e^s ) ,  where  F(s)  is  the  transform  of  f(t)  in  Equation  (14b).  Notice  from 
Equation  (lid)  that  the  function  P(s)  also  contains  the  term  (1  -  e^s )  so  that  in  Equation  (18)  this  term 
can  be  canceled  out  of  numerator  and  denominator.  For  evaluating  Equations  (23),  Qj(s)  and  Q2(s)  may  be 
obtained  by  taking  the  Laplace  transform  of  Equation  (14b),  substituting  s  =  x  +  iy  ,  and  separating  into 
real  and  imaginary  parts  to  obtain  QiR(x,y),  Qu(x,y)  ,Q2R(x,y)  and  Q2i(x,y).  The  resulting  expressions 
are  rather  long  and  are  omitted  here.  The  functions  PR(x,y)  and  Pj(x,y)  are  likewise  obtained  from 
Equation  (1  Id).  The  software  MACSYMA  may  be  useful  in  performing  the  complex  algebra. 

Indicial  Response  Results:  Analytic  Form 

The  indicial  response  results  are  shown  in  Figure  8a  using  the  method  of  Reference  [11]  along  with 
Wagner's  function.  In  the  caption,  a  is  a  variable  used  to  calculate  the  real  value,  a,  of  the  Laplace 
variable  of  integration,  N  is  the  number  of  terms  in  the  Fourier  series,  tmax  is  the  largest  time  for  which 
the  indicial  response  is  the  approximated,  and  E  is  the  relative  error  also  used  to  calculate  the  parameter,  a. 
Each  of  these  quantities  is  a  required  input.  It  has  been  found  that  these  parameters  can  be  varied 
significantly  and  still  provide  good  agreement  with  the  theoretical  value.  For  example  setting  NT=100, 
tmaX~50  ,  and  E=.0001  yields  virtually  the  same  results.  It  is  of  course  necessary  to  choose  a  to  be  larger 
than  the  largest  real  pole  of  y  (s).  The  indicial  response  is  underpredicted  at  the  point  of  discontinuity  at  t=0 
by  one -half  the  actual  value  which,  as  has  been  mentioned,  is  implicit  to  the  method.  For  times  above 
about  one  semichord  the  agreement  between  the  numerical  results  and  the  Wagner  function  is  excellent. 
The  oscillations  in  the  numerical  results  arise  from  the  harmonic  series  used  in  the  approximation.  Shown 
in  Figure  8b  are  the  results  using  the  inversion  method  of  Reference  [12].  Shown  in  the  caption  is  the 
required  input  to  the  algorithm.  The  method  gives  highly  oscillatory  results  at  small  times.  In  fact,  the 
data  below  about  0.2  semichords  have  been  clipped  because  the  amplitude  of  oscillation  is  off  the  scale  of 
the  Figure.  For  larger  times  the  agreement  is  very  good.  These  and  subsequent  calculations  have  been 
done  on  an  IBMp.c.  (Pentium)  using  BASIC  and  requires  approximately  5  minutes  of  computation  time. 
The  calculations  of  the  results  of  Figures  8a  and  8b  require  approximately  thirty  seconds. 

Transform  of  Generalized  Coordinates:  Discrete  Form 

The  treatment  of  q^(t),  q2(t),  and  p(t)  in  the  form  of  discrete  data  points  is  directed  towards  the  situation 
wherein  these  values  would  result  from  experimental  measurements  as  described  below.  The  functions  for 


Time  (semichords) 


Figure  8a  and  8b.  Indicial  Response  Results  for  Analytical  Input  using  8a)  Ref.  [11] 
a=0.01,  N=200,  E=le-6,  trnax=100,  8b)Ref.[12],  a=0.01,  N=200,  ^rnax- 
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qj(t)  and  Q2(0  given  by  Equation  (14b)  and  p(t)  given  by  Equation  (lib)  have  been  used  to  generate  discrete 
values  for  these  functions  at  a  time  interval  of  At  =  0.01  semi  chords  (.0012  sec).  This  time  interval  was 
chosen  because  it  is  representative  of  the  Ohio  U.  tow  tank.  These  data  when  plotted  would  then  closely 
resemble  the  curves  in  Figures  4, 5c,  and  5e.  The  Laplace  transform  of  any  function  f(t)  is  defined  by: 


F(s)  = 


e~stf(t)dt 


(27) 


If  f(t)  is  in  the  form  of  N  discrete  data  points  fQ  ,  f^  ,  ....  at  times  tQ  ,  t^  ,  ....  t^  ,  respectively,  a 
trapezoidal  rule  approximation  to  (27)  is: 


F(s)  =  At  [  i-(  f0  +  e'stN  fN  )  +  J1  e'stj  f  j]  (28) 

2  j=l 

Substituting  s  =  x  +  iy  and  separating  into  real  and  imaginary  parts  yields  F(x,y)  =  FR(x,y)  +  iFj(x,y), 
where: 


FE(x,y)  =  4t  [  i-(  f0  +  e"x,N  cos(jIn)  fN  )  +  2  e"”j tos  Wj)  f  j] 
2  j=l 

N-l 

Fi(x,y)  =  At  [ -  i-e'xtN  sin(ytK)  fN )  -  £  extj  sin  (ytj) f  ■] 

2  j=l 


(28a) 


where  the  identity  e”'-1  =  cos(yt)  -  i  sin(yt)  has  been  used.  For  evaluating  Equations  (28a)  using  the 
method  of  Reference  [11]  (Equation  (22))one  w  ould  set  x  =  a,  and  y  =  toct/T,  whereas  using  the  method  of 

0. 

Reference  [12]  (Equation  (26c))  x  =  a,  and  y  =  -1-cot—  .  It  has  been  found  in  the  present  case  that 

applying  Equations  (28a)  over  a  time  interval  of  approximately  twenty  semichords  gives  an  adequate 
approximation  to  the  Laplace  transforms  of  Qir(x,)0,  Qi j(x,y)  »Q2R(x»y)»  Q2l(x*y)’  ^R^y)  Pi(xJ)- 


Indicial  Response  Results:  Discrete  Form 

The  inverse  results  using  method  of  Reference  [11]  is  shown  in  Figures  9a  and  9b,  where  it  is  seen 
that  the  agreement  with  the  Wagner  function  is  very  good.  A  comparison  of  Figure  9a  and  9b  shows  the 
effect  of  the  number  of  terms  used  in  the  Fourier  series,  with  the  increase  in  the  number  of  terms  increasing 
the  frequency  content  of  the  results.  Figure  9c  shows  the  results  using  the  inversion  routine  of  Reference 
[12]  where  again  there  are  significant  oscillations  near  the  origin. 


Inverse  Formulation  in  Terms  of  Experimental  Data 

Typical  experimental  data  may  include  force  measurements  and/or  accelerometer  measurements.  The 
force  measurements  could  be  either  strain  gauge  or  integrated  pressure  measurements.  Below  we  assume 
the  existence  of  a  set  of  discrete  strain  gauge  normal  force  data  C^t)  and  a  concurrent  set  of  chord  normal 
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Figure  9a  and  9b.  Indicial  Response  Results  for  Discrete  qi(t),  q2(t),  and  p(t)  Input  using 
9a)  Ref.[ll],  cx=0.01,  N=100,  E=le-6,  t 


;max=100,  9b)Ref.[ll],  a=0.01,  N=300, 
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accelerometer  data  a(t)  at  some  point  on  the  structure.  To  initially  investigate  the  use  of  experimental  data, 
while  also  making  use  of  the  analytical  information  above  for  guidance,  we  generate  fictitious  normal  force 
data  and  accelerometer  data.  The  strain  gauge  data  are  the  sensible  force  data  shown  previously  Figure  6c 
(labeled  "Analytical")  and  given  by  Equations  (15a)  and  (15b).  The  accelerometer  data  are  shown  in  Figure 
1 0  and  have  been  generated  using: 

a(t)  =  Diq1(t)  +  D2q2(t)  (29a) 

wbereDj  andD2  are  simply  the  values  of  the  first  and  second  mode  shapes,  respectively,  at  the  point  on 
the  structure  at  which  the  acceleration  is  measured,  and  the  derivatives  of  q^(t)  and  q2(t)  are  obtained  by 
differentiating  Equation  (14b).  The  acceleration  data  is  chosen  to  be  the  chord  normal  acceleration  of 
element  13  measured  on  the  pitch  axis  (see  Figure  lb).  The  force  and  the  accelerometer  data  have  been 
generated  at  a  time  interval  between  data  points  of  0.01  semichords  (.0012  sec),  which  is  reasonable  from  an 
experimental  standpoint.  Note  that  the  analytical  displacement  of  any  point  on  the  structure  is  also  known 
from  the  solution  given  by  Equation  (14b). 

The  displacement  of  the  structure,  d(t),  may  theoretically  be  obtained  by  numerically  integrating  the 
accelerometer  data  This  has  been  done  using  both  a  trapezoidal  method  as  well  as  Simpsons  three  point 
method  Both  methods  have  proven  to  give  unsatisfactory  results  in  that  the  error  in  the  integrated  results 
grows  continuously  as  shown  in  Figure  11a  along  with  the  analytical  displacement  for  comparison.  The 
time  step  between  data  points  was  reduced  to  0.005  with  the  same  result.  No  systematic  attempts  have 
been  made  to  pinpoint  the  source  of  the  error. 

The  displacement  may  alternatively  be  computed  from  the  accelerometer  data  by  numerically  solving 
the  second  order  differential  equation: 

d(t)  =  a(t)  (29b) 

subject  to  the  boundary  conditions  d(0)  =  0,  d(tM)  =  dgs,  where  tM  is  some  time  at  which  steady  state  is 
assumed  to  exist,  and  dss  is  the  steady  state  displacement.  Notice  that  the  steady  state  displacement  can 
be  computed  from  the  steady  state  values  of  q^(t)  and  q2(t)  and  their  respective  mode  shapes  without  regard 
to  actually  solving  Equations  (3).  The  steady  value  of  q^(t)  can  be  computed  directly  from  Equation  (3a) 
and  Equation  (5).  The  only  term  on  the  LHS  of  Equation  (3a)  which  survives  in  the  steady  state  is  the 
stiffness  term,  whereas  on  the  RHS  the  only  surviving  term  is  the  leading  term  in  the  integral  of  Equation 
(5)  which  becomes  CNaLj  ps  s ,  where  ps  s  is  the  steady  value  of  p(x).  The  steady  state  value  of  q^(t)  may 
be  easily  found  from  these  two  surviving  terms.  The  steady  value  of  q2(t)  may  likewise  be  found  from 
Equation  (3b).  Thus,  it  is  reasonable  to  assume  that  ds  s  may  be  known  a  priori. 

Equation  (29b)  has  been  solved  using  a  fourth  order  Runge-Kutta  method  [10]  and  the  results  are  shown 
in  Figure  1  lb  with  the  analytical  displacement  for  comparison.  The  slight  discrepancy  between  the  two  is 
in  part  due  to  the  assumption  of  steady  state  after  25  semichords  in  the  numerical  results,  while  the 
analytical  steady  state  is  reached  assymtotically.  An  advantage  of  the  Runge-Kutta  method  is  that  it  forces 


d(t)  (semichords) 


Time  (semichords) 


Figures  1  la  and  lib.  Displacement  Results  using  1  la)  Simpsons  Method,  1  lb)  Fourth 
Order  Runge-Kutta  Method. 
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the  experimental  displacement  to  the  proper  steady  state. 

The  normal  force  data,  C^t),  and  displacement  data  are  related  to  the  generalized  coordinates  by: 


CN(t)  =  Kiq1(t)  +  K2q2(t) 

(29c) 

d(t)  =  D1q1(t)  +  D2q2(t) 

(29d) 

where  the  constants  Kj  (=12.014)  and  K2  (=5. 148)  can  be  deduced  from  Equations  (15a)  and  (15b).  Taking 
the  Laplace  transform,  letting  s  =  x  +  iy,  and  solving  for  the  real  and  imaginary  parts  of  Q1R,  Qlt,  Q2R, 
andQ2!  yields: 


Q1R(x,y)  =  ( D2CNR(x,y)  -  K2DR(x,y) )/(  D2Ki  -  )  (30a) 

Qll(x,y)  =  ( D2CNi(x,y)  -  K2E>i(x,y) )/(  D2Ki  -  DiK2  )  (30b) 

Q2R(x,y)  =  ( CNR(x,y)  -  K!Q1R(x,y)  )/K2  (30c) 

02i(x,y)  =  ( C\Tj(x,y)  -  K1Qn(x,y)  )/K2  (30d) 


where  C^(x9y)  and  CN|(x,y)  are  the  real  and  imaginary  parts,  respectively,  of  the  transformed  normal  force 
data,  and  DR(x,y)  and  Dj(x,y)  the  real  and  imaginary  parts  of  the  displacement  data,  both  pairs  being 
evaluated  numerically  using  Equations  (28a).  Equations  (30)  provide  the  relations  for  the  generalized 
coordinates  necessary  to  evaluate  Equations  (23).  Figure  12  is  a  flow  chart  of  the  procedure  for  using  the 
experimental  data  in  the  inversion  routine  of  Reference  [11]. 

It  should  be  noted  that  theoretically  it  is  possible  to  arrive  at  the  transforms  of  the  generalized 
coordinates  directly  from  the  accelerometer  data  without  solving  Equation  (29b)  to  obtain  displacement. 
The  Laplace  transform  of  Equation  (29a)  is: 

A(s)  =  s2(  DxQi(s)  +  D2Q2(s)  )  (31) 

where  A(s)  is  the  Laplace  transform  of  a(t).  Setting  s  =  x  +  iy,  Equation  (31)  may  be  combined  with 
Equation  (29c)  to  provide  alternative  relations  for  Q1R>  Q^,  Q2R,  and  Q2I  in  terms  of  CjvjR(x,y), 
CM(x,y),  and  AR(x,y),  Aj(x,y)  where,  as  before,  these  quantities  are  evaluated  using  Equations  (28a). 
This  method  has  proven  unsuccessful  since  Equation  (28a)  is  a  trapezoidal  rule  integration  and  as  has  been 
seen  (Figure  1  la),  integration  of  the  acceleration  data  by  this  method  is  unstable. 

Inverse  Indicial  Response  Results:  Experimental  Data  Formulation 

The  inverse  results  using  the  fictitious  strain  gauge  and  displacement  data  are  shown  in  Figure  13. 
These  results  are  based  on  the  method  of  Reference  [1 1].  The  numerical  results  are  slighdy  larger  than  the 
theoretical  indicial  response.  This  difference  may  be  attributed  to  error  introduced  in  the  calculation  of  the 


I.  Input:  Necessary  Constants 
Input:  Cj^j(t),  d(t),  and  p(t)  versus  time  data 
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II.  Set:  a=.01 
E=lE-6 


N=100 
T  =  0.8  tjjjax 
a  =  a  -  ln(E)/2T 


III.  Set:  x  =  a 

yk  =  krc/T  ,  k  =  0,  1 . N 


IV.  For  each  value  of  (x,yk) ,  k  =  0,  1....N  perform  the  following: 

1.  Compute  B  1R(x,yk),  B2R(x,yk), . B6R(x,yk)  and 

Bll(x,yk),  B2I(x,yk), . Bgjfx.y^  using  Equations  (19) 

2.  Compute  CyR(x,yk),  Cj<r[(x,yk),  DR(x,yk),  Dj(x,yk),  i’R(x,yk),  and 
Pl(x,yk)  using  Equations  (28a) 

3.  Compute  QiR(x,yk),  Qn(x,yk),  Q2R(x,yk),  and  Q2I(x,yk) 
using  Equations  (30) 

4.  Compute  y  R(x,yk)  and  y  i(x,yk)  using  Equations  (23) 


V.  For  Each  Time,  t,  of  Interest  Compute  T(t)  using 
Equation  (22) 


Figure  12.  Flow  Chart  for  Computing  Inverse  Laplace  Transform  using  Strain  Gauge  and 
Displacement  Data. 
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Figure  13.  Indicial  Response  Results  for  Strain  Gauge  and  Displacement  Data  Input  using 
Ref. [11],  a=0.01,  N=200,  E=le-6,  tmax=100. 
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displacement  data  from  the  accelerometer  data  as  has  been  described  in  relation  to  the  results  of  Figure  lib. 
Nevertheless,  the  agreement  is  still  very  good. 

Part  III:  Experimental  Studies 

This  section  incorporates  the  use  of  experimental  strain  gauge  data  in  the  inverse  aeroelastic  model 
described  above.  The  strain  gauge  data  were  taken  as  described  in  the  introduction  to  this  report  and  in 
detail  in  Reference^]. 


Simplified  Aeroelastic  Model 

The  aeroelastic  model  described  in  Parts  I  and  II  of  this  report  has  been  simplified  by  considering  the 
NDOF1  mode  only.  This  simplifies  the  formulation  in  that  the  generalized  coordinates  for  the  NDOF2 
and  TDOF  are  no  longer  required.  This  yields  the  Laplace  transformed  aeroelastic  reciprocal  equations  for 
the  indicial  response  y  (s)  and  the  NDOF1  generalized  coordinate  Qj(s): 


and 


r  x  Bi(s)Qi(s)  +  B2(s)P(s) 
^  ^  B3(s)Qi(s)  +  B4(s)P(s) 


Qi(s) 


[Al(s)  +  A2(s)  y(s)]  P(s) 
A3(s)  +  A4(s)y(s) 


(32) 


(33) 


where  A-^(s)-A^(s)  and  Bj(s)-B4(s)  are  functions  of  the  Laplace  variable  only  and  may  be  deduced  from 
Equations  16  and  17  or  alternatively  from  the  Equations  3  through  10  by  neglecting  the  NDOF2  and  the 
TDOF  coordinates.  By  separating  the  above  Equations  into  real  and  imaginary  parts  as  has  been 
described,  either  Equation  may  be  inverted  to  the  time  domain  to  obtain  either  T(t)  or  q^(t)  using  the 
numerical  procedure  outlined  in  the  previous  section.  For  introducing  the  experimental  strain  gauge  data 
in  Equation  32,  q^(t)  is  obtained  using  Equation  29c  (neglecting  and  p(t)  is  the  experimental  angle  of 
attack.  The  simplified  model  eliminates  the  need  for  a  second  measurement  for  obtaining  q2(0  as  in 
previous  section.  In  Equations  32  and  33,  the  Laplace  transforms  Q^(s),  y  (s),  and  P(s)  are  computed 
numerically  using  Equations  28. 

Figure  14  illustrates  the  effect  of  neglecting  higher  modes  on  the  sensible  force  at  the  strain  gauge. 
The  solid  curve  illustrates  the  sensible  force  considering  higher  modes  as  heretofore  described  and  the  broken 
line  illustrates  the  simplified  model.  As  expected,  the  simplified  model  does  not  exhibit  higher  frequency 
coupling.  The  two  curves  would  be  expected  to  have  similar  (numerical)  Laplace  transforms. 


Validation  of  Reciprocal  Model 

In  this  section  experimental  strain  gauge  data  are  used  to  compute  Qj(s)  and  angle  of  attack  data  are 
used  to  compute  P(s)  to  obtain  the  indicial  response  T(t)  via  Equation  32  and  the  numerical  inversion 
method  of  Ref erence[ll].  This  indicial  response  solution  is  used  to  compute  y  (s),  which  is  then  used  in 
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Equation  33  to  reproduce  the  reciprocal  solution  for  q^t).  This  reciprocity  provides  a  validation  of  the 
model  and  a  measure  of  aeroelastic  effects. 

Figure  15a  illustrates  angle  of  attack  data  for  a  step  onset  of  2°.  The  initial  onset  angle  has  been 
subtracted  off  so  that  the  initial  angle  is  shown  as  zero.  The  Figure  illustrates  two  separate  runs  with  one 
step  up  (increasing  a)  and  one  step  down.  Here  attention  is  focused  on  the  step  up  case.  Figure  15b 
illustrates  qj(t)  data  computed  from  the  strain  gauge  data  using  Equation  29c.  The  initial  values  have 
again  beep  subtracted  off.  Figure  16a  illustrates  the  numerical  solution  for  the  indidal  response  as  well  as 
the  indicial  response  computed  directly  from  the  strain  gauge  data  using:  T(t)  =  C^(t)/Aa  .  Here  the 
numerical  indicial  response  is  multiplied  by  the  steady  state  slope  C^a  (=  8.0/rad)  which  was  measured 
independently  in  steady  alpha  tests.  It  can  be  seen  that  after  about  one  semichord  the  numerical  solution 
tracks  the  mean  of  the  experimental  results.  The  experimental  result  contains  some  higher  frequency 
oscillations  which  are  not  completely  understood.  The  frequency  of  the  oscillations  is  not  close  enough  to 
the  computed  NDOF2  frequency  to  lead  one  to  the  conclusion  that  these  are  associated  with  the  NDOF2 
bending  mode.  The  oscillations  may  be  associated  with  roughness  and/or  misalignment  in  the  tow  tank 
track.  Nevertheless,  the  results  indicate  the  the  effects  of  aeroelasticity  on  the  strain  gauge  output  are 
predominately  confined  to  a  small  period  shortly  after  the  step  input.  Figure  16b  illustrates  the  reciprocal 
numerical  result  for  q^(t)  obtained  using  the  numerical  indicial  response  solution  of  16a.  The  reciprocal 
solution  is  in  agreement  with  the  experimental  input.  The  numerical  result  does  not  contain  some  of  the 
higher  frequency  content  present  in  the  experimental  data  which  is  most  likely  due  to  the  limited  number  of 
terms  used  in  the  Fourier  inversion  routine.  It  is  worth  noting  in  Figure  16a  that  the  numerical  inversion 
routine  predicts  half  the  initial  value  at  time  zero,  which  is  a  consequence  of  the  formulation  [11]. 

Figures  17a  and  17b  illustrate  the  reciprocal  pair  using  the  experimental  data  of  Figures  15a  and  15b 
for  the  step  down  case.  As  seen  in  Figure  17a,  the  effects  of  aeroelasticity  are  confined  to  times  below  five 
semi  chords. 

Central  Difference  Formulation  for  the  Indicial  Response 

Ideally  the  indicial  response  for  the  step  up  and  step  down  cases  would  be  identical.  Comparisons  of 
Figures  16a  and  17a  indicates  some  differences  between  the  two  cases.  These  may  be  due  in  part  to 
slightly  different  onset  conditions  and  aerodynamic  reactions  to  the  step.  In  each  case  the  step  height  was 
approximately  one  degree  and  the  nondimensional  pitch  rate  was  about  0. 15.  The  effective  angle  of  attack 
at  the  leading  edge  in  the  step  up  and  step  down  cases  during  the  pitch  motion  are  significantly  different  due 
to  the  wash  induced  by  the  pitch  motion.  Subsequent  boundary  layer  and  wake  formation  may  be 
influenced  by  these  events. 

To  account  for  differences  in  step  up  and  down  results,  the  indicial  response  has  been  computed  using  a 
central  difference  given  by: 


F(t)  —  ru(t)Actd  +  F(j(t)Aau 


(34) 


Aau  +  Aad 


Figure  15a.  Angle  of  Attack  Data  for  Step  Up  and  Step  Down  Runs  for  an  Onset  Angle  of 
2°. 


Figure  15b.  NDOF1  Results  from  Strain  Gauge  Data  for  Step  Up  and  Step  Down  Runs  for 
an  Onset  Angle  of  2°. 


Figure  16a.  Numerical  Indicial  Response  Results  and  Directly  Computed  Indicial  Response 
for  Step  Up  at  an  Onset  Angle  of  2°. 
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Figure  16b.  Reciprocal  NDOF1  Generalized  Coordinate  Results  and  Directly  Measured 
NDOF1  Data  for  Step  Up  at  an  Onset  Angle  of  2°. 


Time  (semichords)  • 

Figure  17a.  Numerical  Indicia!  Response  Results  and  Directly  Computed  Indicial  Response 
for  Step  Down  at  an  Onset  Angle  of  2°. 


Time  (semichords) 


Figure  17b.  Reciprocal  NDOF1  Generalized  Coordinate  Results  and  Directly  Measured 
NDOF1  Data  for  Step  Down  at  an  Onset  Angle  of  2°. 
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where  the  subscript  "uH  refers  to  step  up  and  "d”  to  step  down.  The  central  difference  results  are  shown  in 
Figure  18  along  with  the  theoretical  result  of  Wagner.  The  present  result  indicates  a  smaller  indicial 
response  than  that  predicted  by  Wagner. 

Figure  18  indicates  that  the  effects  of  aeroelasticity  are  strongest  below  about  two  semi  chords  after  the 
step.  This  is  due  to  the  large  inertial  loads  generated  by  the  rapid  step  input.  This  is  an  important  result 
when  considering  that  it  is  not  necessary  to  experimentally  impart  rapid  step-like  input  motions  as  has  been 
done  in  the  present  tow  tank  studies  in  an  attempt  to  directly  measure  the  indicial  response.  It  should  only 
be  necessary  to  consider  small  amplitude  motions  which  are  locally  linear  in  the  indicial  response  (i.e.  the 
response  is  the  same  over  the  motion  amplitude).  Such  a  motion  might  include  a  more  gradual  small 
amplitude  ramp.  This  type  of  motion  would  not  involve  the  large  inertial  loading  experienced  in  the 
present  rapid  step  motions  and  may  give  better  results.  The  use  of  a  rapid  step  motion  was  motivated  in 
earlier  tests  designed  to  directly  measure  the  indicial  response  which  is  defined  for  a  step  input  The  present 
aeroelasdc  formulation  does  not  require  a  rapid  step  input.  The  present  method  for  determining  indicial 
responses  would  appear  to  be  applicable  to  wind  tunnel  rigs  and  three  dimensional  rigs  so  long  as  the 
aeroelasdc  equations  are  valid. 

High  Angle  of  Attack  Case 

This  section  describes  the  treatment  of  experimental  step  up  and  step  down  cases  for  an  onset  angle  near 
17°.  The  principal  difficulty  in  the  higher  angle  case  is  the  validity  of  the  aerodynamic  loading  terms  in 
the  equations  of  aeroelasticity,  and  in  particular  the  so-called  noncirculator}7  terms.  The  discussion  below 
deals  with  a  simple  method  for  estimating  the  noncirculatory  loads  using  one  of  the  modal  aeroelasdc 
equations.  Below  attention  is  focused  on  the  NDOF2  Equation  3b  which  is  of  the  form: 

M2H2  +  K2q2  =  cxp  +  c2p  -  c^q\  -  c4q2  +  circulatory  terms  (35) 

where  c^  through  c4  are  constants  and  the  circulatory  terms  are  formulated  in  terms  of  the  convolution 
integral.  In  the  linear  case  the  constants  cx  through  c4  all  contain  the  two  dimensional  noncirculatory 
aerodynamic  loading  term,  multiplied  by  structural  terms  (mode  shapes)  which  arise  in  formulating  the 
uncoupled  generalized  Equation  (see  Equations  5-7).  If  one  of  the  constants  is  known,  each  of  the  other 
constants  may  be  computed  therefrom  from  the  known  structural  properties. 

The  constant  c4  may  be  combined  with  M2  and  thus  the  term  apparent  mass.  The  free  vibration 
response  is  therefore  governed  by: 

(M2  +  c4)q2  +  K2q2  =  0  (36) 

which  has  a  harmonic  solution  with  a  frequency  of  vibration  given  by: 


(37) 


Time  (semichords) 


Figure  18.  Numerical  Indicial  Response  Results  for  a  Step  Onset  Angle  of  2  Using 
Central  Difference  Formulation. 
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Thus,  assuming  the  structure  is  linear  so  that  K2  and  M2  are  constants,  a  measurement  of  a>2  would  allow 
one  to  compute  the  the  apparent  mass  c^  using  Equation  (37).  The  remaining  constants  c^  through  C3  may 
be  computed  from  c^.  The  algebra  is  left  to  the  interested  reader. 

An  accelerometer  located  on  mass  element  13  has  been  used  to  measure  the  frequency  response  of  the 
test  rig  at  various  step  onset  angles  of  attack.  A  FFT  was  performed  on  the  accelerometer  signal  and  the 
power  spectral  density  has  been  computed.  The  power  spectral  density  is  simply  the  magnitude  of  the 
complex  Fourier  coefficients  [7].  The  results  for  onset  angles  of  2° ,  1 1°,  and  25°  are  illustrated  in  Figure 
19.  Recalling  that  the  NDOF2  aeroelastic  frequency  based  on  linear  airfoil  theory  was  around  12 
rad/ semichord,  it  is  seen  that  the  NDOF2  response  occurs  at  nearly  the  same  frequency  in  all  three  cases. 
This  suggests  that  the  apparent  mass  and/or  noncirculatory  loading  at  higher  angles  of  attack  may  be  similar 
to  that  for  the  linear  low  angle  case.  These  results  are  preliminary  and  presented  here  as  conjecture. 

Using  the  same  governing  equations  as  in  the  linear  case,  and  analyzing  experimental  data  for  an  onset 
angle  of  17°  in  a  way  parallel  with  that  described  above  for  the  2°  onset  case  yields  the  indicial  response 
shown  in  Figure  20.  These  results  are  based  on  the  central  difference  formulation  above  involving  both  a 
step  up  and  a  step  down  run.  Again  it  appears  that  aeroelastic  reactions  ar  strongest  for  times  below  about 
two  semichords. 


8.00  -i 


51 


6.00  - 


Power 

Spectral 

Density 


4.00  H 


2.00 


NDOF2 


0.00  i — i — i — I — r-ph — I — I — I — I — r— r — i — i — | — i — r  i  i 

0.00  10.00  20.00  30.00 


40.00  '  50.00 


Nondimensional  Frequency  (rad/semichord) 


Figure  19.  Power  Spectral  Density  Distribution  for  High  and  Low  Onset  Angles. 


Figure  20.  Numerical  Indicial  Response  Results  for  a  Step  Onset  Angle  of  17°  Using  a 
Central  Difference  Formulation. 
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Conclusion 

The  research  described  herein  has  focused  on  two  areas.  The  first  is  a  detailed  aeroelastic  analysis  of  an 
airfoil  test  rig  undergoing  rapid  motion  at  low  incidence.  The  analysis  has  led  to  a  relatively  general 
method  to  solve  the  differential  equations  of  aeroelasticity  in  closed  form  using  Laplace  transforms  and 
state-of-the-art  computer  software(MACSYMA)  for  performing  the  required  algebraic  operations.  The 
results  indicate  that  aeroelastic  reactions  may  have  a  significant  effect  of  strain  gauge  data  and  should  be 
considered  in  experimental  design  of  the  type  described  as  well  as  subsequent  data  reduction. 

The  second  area  of  research  has  been  directed  towards  developing  a  method  to  determine  the  airfoil 
indicial  response  from  typical  aerodynamic  measurements  such  as  strain  gauge  and  accelerometer  data. 
The  method  is  essentially  a  solution  of  the  equations  of  aeroelasticity  expressed  in  reciprocal  form  for  the 
indicial  response  given  the  aeroelastic  response.  Because  the  method  is  directed  to  the  reduction  of 
experimental  data  in  discrete  form,  a  numerical  method  has  been  used.  The  method  provides  a  theoretically 
sound  approach  for  determining  an  indicial  response  experimentally,  even  though  the  indicial  response  is 
defined  for  an  instantaneous  step  of  diminishing  amplitude.  It  has  been  shown  that  given  an  aeroelastic 
response  computed  using  Wagner's  function,  the  inverse  method  correctly  determines  the  indicial  response 
to  be  Wagner's  function.  The  method  has  been  applied  to  a  limited  number  of  experimental  cases  as  well. 
The  method  may  have  applications  for  other  types  of  systems  such  as  heat  transfer  and  pure  structural 
modeling  as  the  governing  equations  are  similar  in  form. 

The  author  believes  that  the  type  of  research  described  in  this  report  can  provide  guidance  in  the  design 
of  future  experiments.  It  is  anticipated  that  an  experimental  apparatus  could  be  designed  to  circumvent 
some  of  the  difficulties  encountered  using  the  present  rig.  For  example,  it  should  be  possible  to  design  a 
structure  for  which  it  is  necessary  to  consider  only  a  single  mode  of  vibration,  as  opposed  to  the  system  of 
three  modes  given  by  Equations  (3a)-(3c).  Also,  the  fact  the  present  airfoil  w?as  pitched  about  the  quarter 
chord  simplifies  the  formulation  for  the  aerodynamic  moments.  Furthermore,  it  was  necessary  to  consider 
the  normal/torsional  coupling  in  the  aeroelastic  analysis  primarily  because  the  center  of  mass  of  the  airfoil 
was  not  on  the  pitch  axis.  It  should  be  possible  to  judiciously  add  mass  to  the  test  rig  so  that  the  center 
of  gravity  would  coincide  with  the  pitch  axis  to  minimize  these  effects.  Finally,  as  has  been  noted,  it  is 
not  necessary  to  experimentally  impart  rapid  step-like  input  motions  as  has  been  done  in  the  present  tow 
tank  studies.  The  rapid  step-like  motions  were  motivated  by  attempts  to  directly  measure  the  indicial 
response  which  is  defined  theoretically  for  a  step  motion.  It  is  only  necessary  to  consider  small  amplitude 
motions  which  are  locally  linear  in  the  indicial  response.  Such  a  motion  might  include  a  very  gradual 
small  amplitude  ramp.  This  type  of  motion  will  not  involve  the  large  inertial  loading  experienced  in  the 
present  rapid  step  motions.  The  present  method  for  determining  indicial  responses  would  appear  to  be 
applicable  to  ,  tow  tunnel  rigs,  wind  tunnel  rigs,  and  three  dimensional  rigs  as  well  so  long  as  the 
aeroelastic  equations  can  be  accurately  formulated. 

Future  work  may  focus  on  using  the  present  approach  to  determining  the  indicial  response  over  a  range 
of  onset  angles  of  attack,  and  using  the  responses  in  a  numerical  convolution  integral  to  predict  large 
amplitude  motions  such  as  ramp  motions  for  which  baseline  experimental  data  exist  for  comparison.  This 
should  further  provide  validation  of  the  indicial  responses. 
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Appendix  A 

This  Appendix  gives  the  values  of  the  constants  appearing  in  Equations  (19a)-(19f). 


bj  =  20.5091 

b9  =  -3.4639 

b17  =  6.7969 

b2  =  113.576 

b1Q  =  -25.4171 

blg  =  15.2134 

b3  =  -16.7068 

bj  j  =  -14.7368 

b19  =  64.7242 

b4  =  -9.6866 

b12  =  49.8633 

b2Q  =  -126.9761 

b5  =3.6241 

b13  =  54.0083 

b21  =  -190.8867 

